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Abstract. We use the numerical renormalization group method (NRG) to investigate 
a single- impurity Anderson model with a coupling of the impurity to a superconducting 
"-^J . host. Analysis of the energy flow shows, in contrast to previous belief, that NRG 

Ch I iterations can be performed up to a large number of sites, corresponding to energy 

differences far below the superconducting gap A. This allows us to calculate the 
impurity spectral function A(u>) very accurately for frequencies \u>\ ~ A, and to resolve, 
in a certain parameter regime, sharp peaks in A(u>) close to the gap edge. 
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^ ! 1. Introduction 

00 

The vast progress in nanofabrication during the last decades made it possible to 
study basic physical effects in a very controlled manner. One example of such highly 
^ ■ controllable devices are quantum dots [1] , which are used, amongst various applications, 
oj ■ for a detailed and very controlled study of the Kondo effect [2, 3, 4, 5], which is one 
of the prime examples of many-body phenomena. Below a critical temperature (Kondo 
temperature Tk) a local moment, provided by the spin of an electron occupying the 
quantum dot, gets screened by reservoir electrons within an energy window Tk around 
the Fermi energy. 

These seminal experimental works on the Kondo effect, together with the possibility 
to engineer reservoir properties, raise the following intriguing question: What interesting 
effects may arise if the local moment in the quantum dot is coupled to superconducting 
leads while parameters like the Kondo temperature or the superconducting gap can 
be adjusted arbitrarily? In a superconductor as described by Bardeen, Cooper and 
Schrieffer (BCS) [6], electrons with opposite spin and momentum form Cooper pairs, 
thereby expelling the lead density of states around the Fermi energy. An energy 
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gap occurs. Obviously, when combining a Kondo quantum dot with superconducting 
reservoirs both effects compete: Screening of the local moment by electrons around the 
Fermi energy against pair formation of the latter. 

This competition has attracted a lot of interest and various methods were used to 
analyze the properties of the system in the above-mentioned limits, see for example the 
references in [7] or [8, 9, 10, 11]. In the early 1990's Satori et al. [12] proposed to apply 
the numerical renormalization group method (NRG) [13, 14] to the problem. Since 
then NRG was used to calculate the ground state and subgap bound state properties 
like position and degeneracy [12, 15, 16] as well as their spectral weight [17, 18]. Also 
dynamic quantities like the spectral function [18] or the Josephson current [19, 20] can 
be calculated with NRG, as was done recently. 

The first main goal of the present paper is to gain insight into the way NRG works 
when applied to a system of one local level coupled to a superconducting lead. NRG is 
a well established method for solving strongly correlated impurity problems. Actually, 
it was invented to solve the Kondo problem [13] and since then has been generalized to 
various schemes involving localized states coupled to fermionic [21] or bosonic [22, 23] 
baths. For a review, see [14]. The key idea of NRG is to discretize the conduction 
band logarithmically. This leads to a chain Hamiltonian with exponentially decreasing 
couplings, the so called Wilson chain. It can be solved iteratively by enlarging the 
system site by site: Due to the decreasing couplings every new site can be treated as a 
perturbation of the old system, thus increasing the resolution with every step. 

It is, however, not at all obvious whether NRG still works for superconducting leads 
at energy resolutions well below the gap. This is because the pairing energy A is the 
same at all energy scales, thus remains a constant on-site contribution also in the above- 
mentioned chain structure. At iterations where the coupling of the Wilson chain reaches 
A, it is not obvious whether added sites still can be understood as a perturbation in 
the iterative process or not, that is whether NRG still works or not. 

In order to address this problem, we analyze in detail the flow of the eigenenergies 
during the NRG procedure. We show that NRG is indeed capable to resolve the 
continuum close to the gap without any restriction on the energy scale of the 
superconducting gap. 

The second goal of this work consists in the calculation of the impurity spectral 
function close to the gap edge at zero temperature. As in the Anderson model (with 
normal leads) a Kondo resonance may form. However, due to the superconducting 
property of the leads a gap opens up around the Fermi energy, cutting the resonance. 
Our main interest lies in the study of the continuum contribution to the spectral 
function close to the gap edge, implying the need of high resolution in that regime. 
Our calculations cover not only the regime A > TV, for which the continuum part of 
the spectral function was studied in [18], but also A <C Tk- In the latter regime we find 
a sharp peak at the gap edge, vastly exceeding the Kondo resonance contribution. We 
expect this to lead to an enhanced linear conductance, as observed in a recent experiment 
[24] with carbon nanotube quantum dots coupled to superconducting leads. They report 
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dramatic enhancement of the linear conductance (only) in the regime A <C Tr. 

We find similar behaviour of the spectral function in the non-interacting case where 
an analytical solution exists. We analyze this solution and compare our NRG results 
against it. We find excellent agreement especially at energies close to the gap, as ex- 
pected from our study of the energy spectrum. 

The paper is organized as follows: In section 2 we introduce the superconducting- 
lead Anderson model (SC-AM) and the NRG method. In section 3 the flow of the energy 
spectrum of the SC-AM is analyzed. NRG is shown to work at resolutions well below 
the energy scale of the superconducting gap. Section 4 discusses the calculation of the 
spectral function. In the last part we conclude and summarize our results. 

2. The model and the method 

In this section we introduce our model for the quantum dot coupled to a left 
and right superconducting reservoir, as well as the method we use, the numerical 
renormalization group method [13, 21]. We perform a Bogoliubov as well as a particle- 
hole transformation. We derive the formulas in a general form. For convenience, the 
discussions in the later sections will be restricted to a real order parameter of the 
superconductors, equivalent to only one lead. In the non-interacting limit the system 
can be understood in terms of a simple single-particle picture. The latter can be solved 
exactly and will serve as a guideline for gaining a deeper understanding of the problem. 

2.1. Superconducting-lead Anderson model (SC-AM) 

To describe a quantum state coupled to two superconducting reservoirs, we consider the 
standard Anderson model (AM) for a local level coupled to two metallic, non-interacting 
reservoirs [25], and add a BCS-type term if a, describing pair formation in the leads. 
This SC-lead Anderson model, to be called SC-AM, is then described by 

H = H dot + H hyb + Hi ead + H/^, (1) 

with 

H do t = ^e d n d(J + Un^n di (2) 

(7 

H hyb = Y Y V A° d ° + 4 c «k ff ) (3) 

l=L,R kcr 

Hlead = Y Y £ kcL C «k ff (4) 
l=L,R kcr 

Ha = - Y E A '( e ^ c JkT c Lk| + e-^Q_ ki Q kT ). (5) 

l=L,R k 

Electrons on the dot with spin a = { j , j } are created by S a and interact via the Coulomb 
repulsion U with each other, n da = d) a d a being the charge operator for spin a. In 
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the hybridization term H hyln the coupling strength Vi between dot and lead states is 
assumed to be real and independent of the wave vector k. c lka creates an electron in 
lead I = L,R, respectively. Hi ea d describes the conduction band of metallic leads. We 
assume an isotropic and linearized dispersion. The density of states is then constant, 
Po = 1/2-D, where the band ranges from —D to D. In the following D = 1 will serve 
as energy unit. The pairing of electrons with opposite spin and momentum (Cooper- 
pairs) is described by Ai is the magnitude (later on often called the gap), <pi the 
phase of the order parameter of the superconductor. For simplicity, we assume left-right 
symmetry, i.e. A/ = A, <pi = ±0/2 and V\ = V, where I — L,R, respectively. Then the 
unitary rotation 

1 ( e - ^/ 4 e^/ 4 
" \ -ze~^ /4 ie"^ /4 

of the reservoir operators yields a real Hamiltonian. 

It is well known [26, 27] that the Hamiltonian of a bulk superconductor [equations 
(4) and (5)] can be diagonalized by a Bogoliubov transformation. Note, though, that 
the SC-AM cannot be understood as an AM with a superconducting lead density of 
states (pa = kfcl/ v 1 e \ + A 2 )- This is because the Bogoliubov transformation explicitly 
depends on k. The hybridization term c\ ka .d a + c^-Qko- would transform to a complicated 
object that cannot be simplified by rotating the d operators of the local dot space. 

2.2. The numerical renormalization group method (NRG) 

In the 1970's, K.G. Wilson came up with a scheme for solving the Kondo problem 
nonperturbatively: the numerical renormalization group (NRG) [13]. Since then it 
was generalized to various schemes, describing localized electronic states coupled to 
fermionic [21] or bosonic [22, 23] baths. The NRG allows thermodynamic and dynamic 
properties of such strongly correlated systems to be calculated at zero as well as at finite 
temperature. We first discuss the method for the AM (A = 0), then for A ^ 0. For 
brevity we apply all NRG transformations to the full SC-AM already in the discussion 
of the AM. 

The key idea of NRG is to discretize the conduction band of the reservoir 
logarithmically. The Hamiltonian can then be transformed to a chain Hamiltonian. 
In this representation, equations (3) to (5) are mapped onto 



cos 4 (fto* d * + h.c.) - sin -(/JoX + h.c. 



(6) 

Hlead = \ (1 + A" 1 ) E E E A ~" /2 & (fLfln+Kr + fl + ljlna) (7) 
l=e,o (J n=0 

Ha = " A E [UULi + h.c.) - Ul.fl, + h.c.)] . (8) 

n 

Electrons on site n in lead I = e,o are created by f nlu . The dot level only couples to 
the zeroth site of the so-called Wilson chain Hi ea d, where the hybridization is given by 
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T = np2V 2 (the factor 2 stems from the two leads). A > 1 is the discretization parameter 
of the conduction band. £ n = (1 - A~ n_1 )(l - A~ 2n ~ r )~ 1/2 (1 - A^ 2 "" 3 )" 1 / 2 « 1 for large 
n. The hopping matrix elements between successive sites of Hi ea d fall off exponentially 
with A~ n / 2 . The resulting energy scale separation ensures that the AM can be solved 
iteratively. The recursion relation reads 

H = 1/VX 



H dot + H hyb - A Yl Sl (/ior/icu + h - c - t 

l=e,o a 

H N+1 = ^XH N + l -{\+ A" 1 ) J2J2t» {flNaflN+la + k.C.) (9) 

l=e,o a 

- A AN/2 E 5> (//iv+ir^Wu + h - c ) , 

l=e,o u 

ere si = ±1 for I — e,o, respectively. The initial Hamiltonian of the system is related 
to the NRG Hamiltonian by H = liniAr^oo A - ^ -1 ^ 2 Hjy. This relation is exact in the 
limit A — > 1 and — > oo. 

Sites are added successively and at each step the enlarged system is diagonalized. 
Each added site then acts as a perturbation of order A^ 1 / 2 on the previous part of the 
chain. Consequently, the typical energy resolution 5 n of the AM at iteration n is given 
by <5« M oc A~ n / 2 . Thus, by choosing the length of the chain large enough (so that 
Ar N l 2 is much smaller than all other energies in the problem), all relevant energy scales 
can be resolved and treated properly. When adding a site to the system, the dimension 
of the Hilbert space gets multiplied by the dimension d of the state space of that site, 
yielding d = 4 for a single fermionic lead (empty, singly occupied (either up or down), 
doubly occupied). Therefore the dimension of the Hilbert space increases exponentially 
with the length of the chain. Wilson proposed a truncation scheme according to which 
only the lowest Ak ep t eigenstates are kept at each iteration, thereby ensuring that the 
dimension of the truncated Hilbert space stays manageable. Recently, it was shown that 
by keeping track of the discarded states a complete, but approximate, basis of states 
can be constructed [28, 29]. This can be used to calculate dynamic properties like the 
spectral function A (see equation (20) below) which rigorously satisfy relevant sum rules 
[30, 31], like / dcuA(uj) = 1. 

Applying the NRG mapping also to the pairing term of the SC-AM (as already 
done above), an on-site contribution appears, see equation (8), constant in magnitude 
for each site. In the limit A~™/ 2 ^> A, this additional term hardly affects the properties 
of the system. But, when A - "/ 2 ~ A, it is not obvious whether the added sites still act 
as a perturbation in the iterative process (9) or not, that is whether the energy scale 
separation still works or not. In section 3 we will show that the separation of energy 
scales does work also at resolutions much smaller than the gap. 
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2. 3. Bogoliubov and particle-hole transformations 

Satori et al. [12] have shown that a computationally more convenient representation of 
the Hamiltonian can be obtained by performing a Bogoliubov- Valatin transformation 
[bin,u = l/v / 2(o' fin,cr + fi n -a)] as wen as a particle-hole transformation [q^o- = 
k,2n,a, -aci,2n-i,-cj = &i,2n-i )£ r]- ^ = -1 represents the dot, thus d a = c_i, CT . Applying 
these transformations to (9), the Hamiltonian reads 




H lead = -(1 + A" 1 ) £ £ A^/ 2 e„ (4 ff Qn + i CT + (12) 

l=eo t7,n=0 

oo 
n=0,cr 

Operators in the new basis will always be denoted by a tilde, e.g. rida = d\d a or 
nino = c| n(T Q ncr . The Q nonconserving property of has been transferred to Hdot- 
This has two useful effects: (i) For the symmetric model not only the ^-component 
of the total spin (per iteration), S z n = |Sin=-i(" , M — ^Zni)> but also the particle 
number Q N = ^ jn= _ 1(T (ni w — 1) is a conserved quantum number (note that this 
definition yields Q = for the Fermi sea together with a singly occupied dot at T = 0). 
Therefore the dimensions of the matrices to be diagonalized at each iteration (and 
therefore the numerical effort) is reduced significantly, (ii) Additionally, in the non- 
interacting symmetric case (U, e d = 0) the Hamiltonian takes a very simple quadratic 
form. We will focus on its exact solution in the next section. In section 3 the resulting 
single-particle picture will serve as a tool to gain a deeper understanding of reasons why 
NRG does work for the SC-AM. 

For simplicity, we use = in the following. Then, the odd channel decouples and 
the problem reduces to an effective one-lead system. The resulting model is equivalent 
to that describing an impurity embedded in a bulk superconductor. 

2.4- Single-particle picture 

Some properties of the system show up already in the non-interacting case, U — 0. 
The Hamiltonian is then of quadratic form and we only have to solve a single-particle 
problem. The NRG Hamiltonian [(10)-(13)] can be diagonalized up to a large number 
of iterations exactly - that is without truncating the Hilbert space. One can use 
the resulting exact solution as benchmark for the NRG result. We obtain very good 
agreement in the energy spectrum, thus confirming that NRG is capable of accurately 
treating superconducting leads. In section 3 the single-particle picture will also serve as 
a tool to gain a deeper understanding of reasons why NRG does work for the SC-AM. 
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Without lack of generality we restrict the discussion to the symmetric case, = 0. 
Then the NRG Hamiltonian only contains quadratic terms of the form ajcv, with 
some fermionic operator. For every iteration N the single-particle Hamiltonian can be 
diagonalized by some unitary transformation T to H N = X^=-i £ j a j a j- Here cxj = Tjfli 
and the eigenstates \rij) = aj|vac) satisfy a Mj\rij) = nj\nj). The many-body eigenstates 
and the energy spectrum follow from the Schrodinger equation 

H\m)=E m \m) t £ m =(V SjuA - E , (14) 

with the many-body eigenstates \m) = \n± . . . njv) and eigenenergies E m which are 
calculated w.r.t. the ground state energy E . In the ground state |0) all single- 
particle levels with energy below the Fermi energy ep = are occupied, thus 
Eq = X]z( £ ;<o) 5/ ' Expectation values of local operators are evaluated easily, e.g. 
(Ola-xaLilO) = Ei U^iOla^m^ = £ %<0) {U^. 

The construction of the many-body spectrum from the single-particle energy levels 
using (14) is illustrated in figure 1 for the SC-AM. Figure 1(a) shows a sketch of a typical 
single-particle spectrum. The single-particle level spectrum consists of a continuum 
above and below the gap, |ej| > A (represented by a discrete set of closely-spaced 
levels), as well as one subgap level with energy < Sq < A, the so-called Andreev 
level. Note that because of the discretized conduction band, we also have a discretized 
continuum. 

The sketch also demonstrates the construction of the lowest lying many-body 
eigenenergies using (14). For T > 0, no single-particle level exists at the Fermi energy 
and the many-body ground state is a singlet (S z = 0, Q = 0). The first excitation 
is a degenerate doublet (£1,2 = £0, S z = ±1/2, Q = 1), corresponding to the bound 
single-particle level, occupied by either a spin up or down. If e < A/2, an additional 
subgap state forms (E 3 = 2sq < A, S z — 0, Q — 2), corresponding to spin up and 
down occupying the subgap single-particle level. Otherwise E 3 > A is part of the 
continuum energies. A concrete example of the many-body as well as the single-particle 
eigenenergies is shown in figure 1(b). Both the single-particle levels (stars, crosses) 
as well as the resulting many-body eigenenergies are shown. On the vertical axis the 
energies of the many-body eigenstates constructed in figure 1(a) are specified. 



3. Energy spectrum 

In this section we analyze the many-body energy spectrum generated during the iterative 
NRG procedure. As already mentioned in the last section, the spectrum consists of a 
continuum above and subgap bound states below the gap A. The competition between 
the Kondo effect and Cooper pair formation is reflected in the ground state properties 
of the system. A detailed analysis of the structure of the continuum with the help of the 
single-particle picture reveals that, interestingly, energy scale separation is even more 
efficient at energy scales smaller than the gap (compared to the AM). 
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Figure 1. Single-particle and many-body eigenenergies for the SC-AM (U = s<i = 0). 
(a) Schematic sketch of a typical single-particle level spectrum of the SC-AM. The 
continuum is represented by a discrete set of closely-spaced levels. The construction 
of the lowest lying many-body eigenstates in terms of single-particle states according 
to (14) is illustrated. The corresponding energies E m (w.r.t. the ground state energy) 
and quantum numbers of the many-body states are also given. Two dots stand for a 
doubly occupied level, (b) Single-particle energies £j (ej > 0: stars, ej < 0: crosses), 
as well as the corresponding many-body eigenenergies E m (lines) for A = 10~ 4 and 
r/A = 0.3 at a late NRG iteration (S n <§C A). The many-body eigenenergies depicted 
in (a) are specified on the right hand side. The single-particle continuum energies 
go like — A oc A 2 -*, see section 3.2. Due to many-particle excitations, this dense 
profile repeats as substructures in the many-body spectrum at mA + reo (m = 1,2,..., 
r = 0,l,2). 



Figure 2 shows the 280 lowest lying many-body eigenstates for the even NRG 
iterations of a SC-AM in different regimes of TV/A. We first discuss the case A = (AM, 
figure 2(a),(d)), then A^0. For the AM the effective level spacing of the Wilson chain 
drops exponentially with every added site (see discussion above). The energy resolution 
of the kept states is enhanced exponentially with increasing iteration n, see figure 2(a). 
Thus, an appropriate way of visualizing the physics at different energy scales is given 
by the rescaled energy spectrum. In these "energy flow diagrams", the eigenenergies 
are plotted in units of A _n / 2 oc 5^ M , see figure 2(d). Only at energy scales where the 
system changes its properties, the flow of the eigenenergies changes. For the AM we are 
interested in the lowest of these scales, the Kondo scale Tk = -y/^f exp [§^§(U + £d)], 
indicated by dashed arrows in the figure. For details of the various fixed points of the 
AM see e.g. [21]. 

In contrast to the exponential decaying couplings of the Wilson chain (7), the 
on-site contribution A of the pairing term (8) is constant in magnitude for each 
site. Consequently, for 5 n < A, the BCS contribution is a relevant perturbation and 
determines the physics of the system. Typical energy spectra (or energy flow diagrams, 
respectively) for finite A are shown in figure 2(b,c) (or 2(e,f)). At energy scale A, the 
exponential reduction of the eigenenergies crosses over to a saturation towards A. The 
characteristic gap as well as the subgap Andreev bound states form. 
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The structure of the continuum energies near the gap will be discussed in section 
3.2. We show there that even though the on-site terms of do not fall off exponentially 
like the couplings of the Wilson chain, the energy scale separation (the heart of NRG) 
still works. 
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Figure 2. Energy spectra (a,b,c) and the corresponding energy flow diagrams (d,e,f) 
for the lowest 280 eigenenergies of the even NRG iterations of a SC-AM. U = 0.2, 
T = 0.01, T K = 1-25 • 10~ 5 , A = 2.5 and A kcpt = 512 in all plots. The iteration 
numbers Nt k or where 5^ M « Tk or A, respectively, are indicated by dashed 
or solid arrows. (a,d): A = 0, AM. Since S^ M oc A - ™/ 2 , the eigenenergies (a) fall 
off exponentially with n and the energy flow diagram (d) converges. At energy scale 
Tk the localized spin gets screened by the conduction electrons and the Kondo singlet 
(ground state) forms. (b,c,e,f): A > 0. At energy scale A the exponential decrease 
crosses over to a saturation towards A. The characteristic gap as well as the Andreev 
bound states form. Consequently, the flow diagram energies grow with AA™/ 2 due 
to rescaling. For comparison, the result for the AM is indicated (brown), too. For 
A <C Tk (b,e) it is energetically favorable to break Cooper pairs and lower the energy 
by Tk by screening the local spin, thus the ground state is a singlet. For A ^> Tk 
(c,f) Cooper pair formation dominates the low energy properties and the ground state 
is a doublet. 



3.1. Competition between Kondo effect and Cooper pair formation 

The system possesses two energy scales determining the characteristics of the system, Tk 
and A. The resulting competition between Kondo effect and Cooper pair formation is 
reflected in the ground state properties of the system: For T K <C A, is the dominant 
term of the Hamiltonian. The system lowers its energy by A by the formation of Cooper 
pairs leading to effective spin zero bosons (singlets) in the reservoir. The lead density 
of states gets depleted, a gap from —A to A forms. Therefore no electrons are available 
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to screen the localized spin. Consequently, the ground state is a spin doublet with 
5*2 = ±1/2. In case when the Kondo effect is dominant (Tk 3> A), it is energetically 
favorable to break Cooper pairs so that the localized spin on the quantum dot gets 
screened by the non-paired electrons near the Fermi energy. The energy is lowered by 
Tk and as ground state the typical Kondo singlet forms (S z = 0). The influences of the 
different scales is also apparent in the NRG energy flow diagrams, see figure 2(e),(f). 
For example, for Tk ^> A, the effect of H A sets in at iterations after the Kondo fixed 
point is reached (i.e. at lower energies). A phase diagram for the singlet and doublet 
ground state including the spectral weight of the bound states was recently derived by 
[18] (also using NRG) for the whole regime of A, V. 



3.2. Analysis of the continuum 

The key feature of NRG is the energy scale separation: The couplings between successive 
sites of the Wilson chain describing a normal lead fall off exponentially, therefore each 
added site can be treated as a perturbation of the previous system. However, when 
generalizing the AM to superconducting leads, a constant on-site energy A is added at 
each site (see (9)). In order to understand why NRG works even at resolutions well 
below A, we now take a closer look at the structure of the continuum produced during 
the iterative NRG procedure. 

We therefore analyze the (positive) continuum of the single-particle problem. Figure 
3 shows an example of a single-particle spectrum for three different scalings of the 
vertical axis: In (a) the (unsealed) spectrum is plotted versus the iteration number. The 
NRG eigenenergies decrease with iteration number n and tend towards A. However, 
the decrease of the continuum eigenenergies depends on whether n is smaller or larger 
than N A , the iteration number for which 5^ M ~ A. We find the following asymptotic 
behaviour: 

/ Ejn oc A j ~ n / 2 for n<N A , 

£jn ' \e' jn = e jn -A oc A 2 ^ n for n > N A . 1 j 

Primed energies will henceforth always be understood to be measured relative to 
A. Both relations of (15) are illustrated in figures 3(b) and 3(c), by plotting the 
eigenenergies Sj n and e'j n in units of A - "/ 2 and A~ n , respectively. 

These results can already be understood by further reducing the problem to only 
the superconducting reservoir (Hi + H A ): For fixed iteration n, the j-dependence of 
(15) reflects the standard logarithmic discretization of the continuous conduction band, 
according to which the single-particle energies of the Wilson chain grow in powers of 
A, i.e. £fc oc A- 7 [21]. Inserting this into the single-particle dispersion relation of BCS 
quasiparticles, the effective discretization of the isolated superconducting reservoir is 
obtained, as sketched in figure 4. The limits yield 

£ k « Ek -> A^ for e k > A (n < N A ), 



^ = V fi 2 + Aa: ± -> A* for £fc «A (n > N A ). (16) 
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The n-dependence of (15) follows heuristically from considering the coupling of two 
neighbouring sites n — 1 and n, with hopping matrix element t n ~ 5^ M and on-site 
energy ±A (see (12) and (13), respectively). Since t n ~ A - ™/ 2 , the eigenstates A± of 
this two-state problem show the asymptotic behaviour 



Note that the A _n//2 versus A~ n scaling of the two limits of (15) implies that 
the energy scale separation is actually more efficient in the second limit, when the 
gap dominates the spectral features. This establishes one of the central results of the 
present paper: The energy scale separation, being the heart of the NRG approach, is 
not impaired but rather enhanced (w.r.t. A) by the presence of the energy gap in the 
superconducting leads. Increasing the chain length leads to an exponential enhancement 
of the resolution at the continuum edge at A. 

The many-body spectrum is constructed according to equation (14). For n < N A , 
as in the AM, the mean level spacing (at fixed iteration) does not depend on energy, 
as can be seen in figure 2(d-f). For n > N A the gap forms, and every many-body state 
with energy E m < A + Eq can only stem from adding one electron (or hole) to the 
Fermi sea. Eq is the energy of the single-particle subgap level. Therefore in this regime 
Ejan = £jn- Due to many-particle excitations, this dense profile repeats as substructures 
in the many-body spectrum at mA + rEo (m = 1, 2, . . .; r — 0, 1, 2), as can be found in 
figure 1(b), which was calculated for some late iteration with n ^> N A . 

4. Spectral function 

In this section we study basic properties of the spectral function that are special for the 
SC-AM. We therefore begin by analyzing the analytic solution of the non-interacting case 
and find not only a continuum for energies \u\ > A together with subgap resonances 
(as expected from the energy spectrum), but also a sharp peak at the gap edge for 
A < T. The agreement between NRG results and the analytic solution is excellent, 
especially when resolving that sharp feature, again confirming that NRG is valid also 
at resolutions well below the gap. Subsequently, we consider spectral functions at finite 
U, which show some of the same feature as found for the non-interacting case. 

Knowledge of the energy spectrum suffices to calculate thermodynamic quantities, 
such as the impurity specific heat. We focus here instead on the more complex 
calculation of the local spectral function A a (oj). As this is a dynamic quantity, all 
energy scales have to be taken into account even at temperature zero. The spectral 
function is defined as 



where G^(u) is the Fourier-transformed of the retarded Green's function G^(t) = 
— i9(t)([d a {i), d£.(0)]+). Motivated by the structure of the eigenspectrum we distinguish 




oc A~ n / 2 for n< N A , 
oc A~ n for n > N A . 



(17) 



AxM = --ImG« 



(18) 
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iteration n 

Figure 3. Positive single-particle energies plotted for three different scalings of the 
vertical axis for even iterations. The vertical dashed line indicates the iteration 
for which 5^ M ps A. (a) Ener gy versus iteration number. The single-particle levels 
flow towards the gap and form a continuum. The bound state can be also seen, (b) 
The eigenenergies e representing the continuum in units of A - "/ 2 , corresponding to 
the customary scaling for energy flow diagrams of the (normal) AM. For n < Na, the 
single-particle energies run horizontally, since they obey Ej n oc A- 7- ™/ 2 , see text and 
[21]. (c) e 1 in units of A - ". The single-particle energies run horizontally for n > Na, 
since here they scale as e'j n oc A 2j ~". The line sloping upwards at the left of figures 
(b) and (c) is due to the fact that at every iteration a degeneracy is split and hence 
an extra eigencnergy is generated. 




Figure 4. Logarithmic discretization of Ek leads to high resolution at the band edge 
of the BCS quasiparticle eigenenergies. 

two contributions to the spectral function, 

A a (uj)= w m 5{u- E m )+A a {uj)6(\uj\- A). (19) 

|£ m |<A 

w m denotes the weight of excitations to the Andreev bound states, which contribute as S- 
peaks within the gap between — A and A. A a (uj) represents the continuum contribution, 
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with \oj\ > A. In the following we present results for the continuum part of the spectral 
function. We focus on the behaviour of the continuum contribution of the spectral 
function close to the gap. 



4.1. NRG 

NRG calculations of the spectral function are based on the Lehmann representation: 

A a (oj) = 1 J2( e ~^ + e "^) I H4K) I 2 S(u - (E m - E m ,)). (20) 

mm' 

Here Z = ^2 m e~ Em ^ kBT is the partition function at temperature T, and \m) denotes an 
exact eigenstate of the Hamiltonian with eigenenergy E m . The matrix elements of these 
operators as well as the eigenenergies can be calculated with NRG at all energy scales. 
The 5-peaks of the continuum contribution are broadened as described by [32, 31]. As 
expected from the findings about the spectra, broadening w.r.t. A (i.e. using id' = \u\ — A 
in [31]) leads to good results, see below. We use the full density matrix NRG [31, 30] 
at effective temperature zero, i.e. at temperature much smaller than any other energy 
scale of the problem. Only matrix elements connecting the ground state(s) with excited 
states then contribute. 

As NRG parameters we choose A = 1.8 and use z- averaging [33] (where for fixed 
A data is averaged for different discretizations) with an interval spacing of Sz = 0.05 
to impoove the results. We keep Ak ep t = 1024 at the first 6 to 10 iterations, and use 
Akept = 512 for the rest of the iterations. We test NRG against the analytic solution at 
U — and find excellent agreement, implying that the broadening procedure as well as 
the choice of Ak ep t are adequate for the present problem. 



4-2. Spectral function for U = 

For analyzing the basic properties of A(u), we first review the non-interacting problem. 
There, the local (retarded) Green's function is known exactly [18]. The continuum 
contribution to G° then reads 

G°aH = {(w + e d ) + i T PA }, (21) 

with D(uS) = {uj 2 — T 2 — e"j) + i 2TupA, and pa the density of states of a bulk 
superconductor. At temperature T = 0, the latter has the limits 

9 (H-A)H f^K), for u/«A, 
\Iuj 2 - A 2 1, for A < 

For A = 0, equation (21) simplifies to the well known formula for the AM, G° a=0 (uj) = 
(id — Ed — ir)^ 1 . For the spectral function we obtain from equations (21) and (18), 

A{U) = (id 2 -e 2 - F 2 ) 2 + (2F W 1 ^ (23) 
This function is shown for various parameter combinations in figure 5. The common 
features are (i) the atomic resonance of width (half width half maximum) ~ T centered 
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e d =0 



e d /r=4 




1e-12 1e-08 0.0001 
co'=co-A 



1e-12 1e-08 0.0001 
co'=co-A 



Figure 5. Continuum contribution A{uS) to the spectral function for U = 0, T = 0.008 
and A = 10~ 4 as well as A = 10~ 2 obtained from equation (23). (a,c) show a 
symmetric and (b,d) an asymmetric SC-AM. In (a,b), a linear scale is used revealing 
sharp peaks near the gap if oj' c <C A. The insets, which zoom in the region of the gap 
edge, show the full height of the near-gap peaks. They also indicate, using position 
and length of arrows, the energy and weight of the subgap contribution for the case 
A = 10 -2 (calculated with NRG). In (c,d), for the same data a log-log scale is used to 
elucidate the asymptotic behaviour of (25) (dashed and dotted lines) and (26) (dashcd- 
dotted lines). 



at Ed, reflecting the level broadening due to the level-lead coupling and (ii) a gap from 
—A to A, with (iii) bound states at some energy ±ujb inside the gap. The energy and 
weight of the subgap contribution is indicated for A = 1CT 2 by position and length of 
the red arrows in the insets of figure 5(a,b). Here they are calculated with NRG, but 
can be also obtained analytically, see equation (7) of [18]. Note that for finite T bound 
states exist also for Ed 3> A, see figure 5(b). They asymptotically approach the gap 
edge for \s d \ — > oo, in accordance with equation (7) of [18]. 

Additionally, the continuum part of the spectral function may feature near-gap 
sharp peaks, point of interest in the following discussion. The behaviour near the gap 
edge can be approximated (using equation (22) and writing s = sign (a;)) by 

A{U) * (A" - P + + 2TW/u> TpA/7T (24) 

(sA 2 + r^ +r2 IW/tt cxv^ 7 , for u/<a4 

(sA+e d ) 2 +F 2 t / 1 r l „ / ^ a \ / 

(Ai-ri-e^+^AZ r >A/ 7r a 7Z7> fOT U c < W < A. 
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co'=co-A co'=co-A 



Figure 6. Continuum contribution A(ui') to the spectral function for [7 = 0, 
r = 0.008, as obtained from (23), for (a) various values of A and = (symmetric 
model) and (b) A = 10~ 4 and various values of Ed (asymmetric model). For uJcAa 
sharp peak forms at the gap edge. The increase (decrease) of A{uj') goes as u/2 (u>'~ 2 ). 
NRG results, shown for A = 10~ 4 only (thick black lines), arc in excellent agreement 
with the exact analytical results for u>l <C A. Note in (b) that the height of the near-gap 
peak at the gap edge does not depend on e d . 



The limits given in (25) are indicated by the thin lines in figure 5(c,d): A(u>) increases 
as Vu? when to' is increased from (dashed lines), decreasing again for to' > uj' c = 
(A 2 -Y^-l^+iY 2 A. 2 ( wmc h is the zero of derivative of equation (24)). If cu' c <C A,T, this 
leads to a very sharp near-gap peak which decreases as pa (dotted lines). Then the near- 
gap spectral function is greatly enhanced compared to the AM (where the symmetric 
case yields A(0)nT = 1). 

The solution of the AM (dash-dotted lines) describes the high-energy limit of the 
SC-AM: 

A(u>) « ^ = --ImG« for A < \lu\. (26) 

(u - e d y + 1 A 71 

The emergence of the near-gap peak is depicted in figure 6(a) for the symmetric 
model, where the gap A is varied over four orders of magnitude, starting from A = 0. 
For A = 10~ 4 , we also show numerical NRG results (fat solid line). Their agreement 
with the analytical results is excellent. The height of the near-gap peak does not depend 
on Ed, see figure 6(b), where Ed is increased up to 8r for fixed A = 10 -4 . 



4-3. Spectral function for finite U 

Figure 7 shows spectral functions of the SC-AM for finite U and —U < Ed < 0. The two 
atomic resonances of width V are now separated by the Coulomb repulsion, thereby they 
are centered near Ed and Ed + U. The Coulomb repulsion also drives the Kondo effect, 
yielding a sharp resonance of width Tk pinned at the Fermi energy. This resonance is cut 
by a gap reaching from —A to A, reflecting the influence of the superconducting lead. 
Depending on the ratio T K / A, the Kondo resonance can be cut completely (T K / A <^ 1) 
or emerge clearly (T K /A 3> 1). 
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Figure 7. Continuum contribution A(ui) to the spectral function for U = 0.6, 
r = 0.049, calculated using NRG. (a,c) show a symmetric and (b,d) an asymmetric 
SC-AM. In (a,b), a linear scale is used revealing the sharp near-gap peaks appearing 
for A <C Tjc- The position Ed of the local level is indicated by an arrow. The insets 
zoom into the gap edge and show the sharp peaks as well as the subgap resonances 
(for A = 3 • 10~ 2 , indicated by arrows). In (c,d), A{u>') is plotted on a log-log scale for 
various parameters. The asymptotic behaviour of the near-gap peaks is in agreement 
with that given in (25) and (26). In particular, for oj' — > 0, the continuum edge 
decreases as \JU . (c) Ed = —U/2, Tk = 10~ 3 , various A. The singlet-doublet 
transition occurs between A/2\- = 0.3 and 3 (orange: ground state is doublet, black: 
singlet), (d) A = 3 • 10~ 5 , various Ed, A <C Tk always. 



Additionally, a similar near-gap feature as found for the non-interacting case may 
emerge. For TV /A ^> 1, a sharp resonance forms at the gap edge, highly exceeding 
the Kondo resonance, see figure 7 for the symmetric case (where the height of the 
Kondo resonance is given by I/71T) as well as the antisymmetric case. The asymptotic 
behaviour of the near-gap peak is in agreement with that given in equations (25) and 
(26). In particular, for cu' — > 0, the continuum edge decreases as y/uS. 

5. Conclusion 

The NRG is a well established method for a variety of quantum impurity models. It 
is usually applicable in the whole parameter range and allows to calculate physical 
quantities for a wide range of temperatures and frequencies. In this paper we showed 
that in the presence of a superconducting reservoir, NRG provides information for 
resolutions far below the energy scale of the gap A. Moreover, Wilsonian energy scale 
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separation, being the heart of the success of the NRG approach, is not impaired but 
rather enhanced by the presence of the energy gap of the superconducting leads. This 
allows sharp features of spectral functions at the continuum gap edge to be resolved. 
Our calculations of the impurity spectral function cover the whole region from A ^> T# 
to A <C Tk- In the latter case, we find a sharp peak at the continuum gap edge, vastly 
exceeding the Kondo resonance contribution. We expect this to result in an enhanced 
linear conductance, as recently reported for experiments with carbon nanotube quantum 
dots coupled to superconducting leads [24]. 

The ability of the NRG to resolve spectral functions at energy resolutions well 
below the gap should be useful for other problems as well. As discussed in detail in [34], 
the problem of an impurity in a superconducting host can be mapped (under certain 
conditions) to a model in which the superconductor couples to a normal metal, with a 
modified density of states. For the problem studied in this paper we have performed 
such a mapping; the resulting Hamiltonian is given in Equation (8). In this case, the 
oscillating on-site energies, (— l) n A, generate a hard gap of width 2A. 

This connection allows us, in principle, to calculate the dynamic quantities for 
impurity models with arbitrary gapped bath spectral function (but for this one still 
has to develop an algorithm which produces directly the chain parameters from the 
hybridization function). Such calculations might also help to improve the resolution 
of the NRG in DMFT calculations for the Hubbard model where the standard 
implementation of the NRG does not describe the shape of the Hubbard bands properly 
(for dynamic DMRG calculations for this problem, see [35]). 
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